{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## with seasonal adjustment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using DataFrames\n",
    "using QuadGK\n",
    "using Distributions\n",
    "using Gadfly\n",
    "using JLD\n",
    "\n",
    "include(\"invwishrnd.jl\");\n",
    "include(\"logsplineconversion_cps.jl\");\n",
    "include(\"vech.jl\")\n",
    "include(\"nearestSPD.jl\")\n",
    "include(\"nearestSPDinv.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set directories\n",
    "datadir      = \"$(pwd())/data/\"\n",
    "crossdensdir = \"$(pwd())/results/\"\n",
    "#crossdensdir = \"$(pwd())/results/no_lbound/\"\n",
    "MECovdir     = \"$(pwd())/results/StateSpace_seasonality/MECov/\";\n",
    "loaddir      = \"$(pwd())/results/StateSpace_seasonality/NoKron/\";\n",
    "savedir      = \"$(pwd())/results/StateSpace_seasonality/NoKron/\";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "114×9 Array{Float64,2}:\n",
       " -3.93868     0.821879  -0.00753001  …  -1.5749    -5.10393    -13.1708  \n",
       " -0.978684    0.461917  -0.00755301     -3.16288   -4.80928     -0.026007\n",
       " -2.24868    -1.64044   -0.00753181     -4.19118   -3.67302      0.340863\n",
       "  1.97132     2.03784   -0.00261831     -3.3511    -2.75116     -3.80237 \n",
       "  0.971316   -1.1044    -0.00752541     -0.667587  -7.30576     -3.68615 \n",
       " -2.45868    -2.71921   -0.00411361  …  -1.9715    -0.163731    -1.42534 \n",
       " -5.12868    -6.24284    0.00117729     -2.68107   -3.42094     -3.23939 \n",
       " -2.79868    -4.51119    0.0126176      -3.0194    -3.14234     -2.19326 \n",
       "  2.44132     0.412119   0.00665049     -0.61075   -7.01013     -6.4938  \n",
       "  0.371316   -0.917889   0.00573779     -1.81668    0.0151863    0.24218 \n",
       " -0.638684   -1.03152    0.0103773   …  -1.22499   -2.34361      6.30853 \n",
       "  8.58132     2.11829    0.0177246      -0.099651  -1.01116      9.03291 \n",
       "  1.07132     1.65404    0.0134918       1.67949   -2.25689      3.21055 \n",
       "  ⋮                                  ⋱                                   \n",
       "  0.0813158  -0.204167  -0.00556071     -1.23002    3.12139    -17.7216  \n",
       "  0.201316    1.14245   -0.00283021      1.8341    -0.674963   -16.0287  \n",
       "  0.141316    0.658817  -0.00686391      0.466523   2.51094    -26.6888  \n",
       " -1.73868    -0.555858  -0.00987281  …   1.97738    1.30416    -21.8155  \n",
       " -3.07868    -1.6616    -0.00999391      1.86887   10.1936      -8.12406 \n",
       " -1.63868    -1.43592   -0.00713391      3.74134    2.92282    -10.3046  \n",
       " -3.07868     0.172584  -0.00965011      2.67052   -2.07881    -24.2432  \n",
       "  2.97132     0.575579  -0.0102608       1.82483    7.20566    -14.5822  \n",
       " -0.808684   -0.402868  -0.012267    …   0.167062   4.11781    -20.4123  \n",
       " -0.938684   -0.774091  -0.0118286       1.88315    8.72911     -2.09123 \n",
       "  0.381316    0.970574  -0.0144698       3.51786    3.3736     -24.4942  \n",
       "  1.75132     0.948956  -0.0165418       1.61441    4.22242    -22.6454  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#include(\"LoadData.jl\")\n",
    "include(\"LoadData_seasonality_select.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "StatePMeanRaw = Array(readtable(savedir * \"StatePMeanData.csv\",header=true))[:,2:end]\n",
    "stateperiod = Array(readtable(savedir * \"StatePMeanData.csv\",header=true))[:,1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "113×9 Array{Float64,2}:\n",
       " -0.978684    0.461917  -0.00755301  …  -2.82722   -3.25569    -0.967999\n",
       " -2.24868    -1.64044   -0.00753182     -3.90041   -2.66708     1.18689 \n",
       "  1.97132     2.03784   -0.00261832     -2.56571   -3.26329     1.12445 \n",
       "  0.971316   -1.1044    -0.00752541     -0.120048  -6.19445    -3.24353 \n",
       " -2.45868    -2.71921   -0.00411361     -2.33852   -1.20084    -2.64953 \n",
       " -5.12868    -6.24284    0.00117728  …  -3.20983   -2.97692    -3.51706 \n",
       " -2.79868    -4.51119    0.0126176      -2.32152   -3.23513     0.999562\n",
       "  2.44132     0.412119   0.00665049     -0.169813  -6.18572    -5.97228 \n",
       "  0.371316   -0.917889   0.00573778     -2.37997   -0.725521   -1.41909 \n",
       " -0.638684   -1.03152    0.0103773      -3.0154    -1.54507     1.124   \n",
       "  8.58132     2.11829    0.0177246   …  -0.934006  -1.05867     6.71228 \n",
       "  1.07132     1.65404    0.0134918       0.544355  -3.14331    -1.60762 \n",
       " -0.298684    0.990175   0.0100448      -2.44934   -1.73303    -3.97468 \n",
       "  ⋮                                  ⋱                                  \n",
       "  0.0813158  -0.204167  -0.00556071     -0.229222   2.74557   -14.1128  \n",
       "  0.201316    1.14245   -0.00283021      1.89075    0.242917  -16.7453  \n",
       "  0.141316    0.658817  -0.00686391      1.86411    1.60736   -21.3516  \n",
       " -1.73868    -0.555858  -0.00987281      1.93302    1.89942   -21.4891  \n",
       " -3.07868    -1.6616    -0.00999391  …   1.11189    7.71084   -10.6303  \n",
       " -1.63868    -1.43592   -0.00713391      2.7218     2.78431   -12.7085  \n",
       " -3.07868     0.172584  -0.00965011      2.67308   -0.620626  -24.0892  \n",
       "  2.97132     0.575579  -0.0102608       1.18247    5.73092   -16.0862  \n",
       " -0.808684   -0.402868  -0.012267        0.361218   4.15119   -18.698   \n",
       " -0.938684   -0.774091  -0.0118286   …   1.31161    7.76325    -4.38414 \n",
       "  0.381316    0.970574  -0.0144698       3.60196    1.29311   -22.8787  \n",
       "  1.75132     0.948956  -0.0165418       1.61441    4.22242   -22.6454  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "statepmean = StatePMeanRaw[:,1:size(data)[2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "113×9 Array{Float64,2}:\n",
       " -0.978684    0.461917  -0.00755301  …  -3.16288   -4.80928     -0.026007\n",
       " -2.24868    -1.64044   -0.00753181     -4.19118   -3.67302      0.340863\n",
       "  1.97132     2.03784   -0.00261831     -3.3511    -2.75116     -3.80237 \n",
       "  0.971316   -1.1044    -0.00752541     -0.667587  -7.30576     -3.68615 \n",
       " -2.45868    -2.71921   -0.00411361     -1.9715    -0.163731    -1.42534 \n",
       " -5.12868    -6.24284    0.00117729  …  -2.68107   -3.42094     -3.23939 \n",
       " -2.79868    -4.51119    0.0126176      -3.0194    -3.14234     -2.19326 \n",
       "  2.44132     0.412119   0.00665049     -0.61075   -7.01013     -6.4938  \n",
       "  0.371316   -0.917889   0.00573779     -1.81668    0.0151863    0.24218 \n",
       " -0.638684   -1.03152    0.0103773      -1.22499   -2.34361      6.30853 \n",
       "  8.58132     2.11829    0.0177246   …  -0.099651  -1.01116      9.03291 \n",
       "  1.07132     1.65404    0.0134918       1.67949   -2.25689      3.21055 \n",
       " -0.298684    0.990175   0.0100448      -2.96456   -2.00878     -5.855   \n",
       "  ⋮                                  ⋱                                   \n",
       "  0.0813158  -0.204167  -0.00556071     -1.23002    3.12139    -17.7216  \n",
       "  0.201316    1.14245   -0.00283021      1.8341    -0.674963   -16.0287  \n",
       "  0.141316    0.658817  -0.00686391      0.466523   2.51094    -26.6888  \n",
       " -1.73868    -0.555858  -0.00987281      1.97738    1.30416    -21.8155  \n",
       " -3.07868    -1.6616    -0.00999391  …   1.86887   10.1936      -8.12406 \n",
       " -1.63868    -1.43592   -0.00713391      3.74134    2.92282    -10.3046  \n",
       " -3.07868     0.172584  -0.00965011      2.67052   -2.07881    -24.2432  \n",
       "  2.97132     0.575579  -0.0102608       1.82483    7.20566    -14.5822  \n",
       " -0.808684   -0.402868  -0.012267        0.167062   4.11781    -20.4123  \n",
       " -0.938684   -0.774091  -0.0118286   …   1.88315    8.72911     -2.09123 \n",
       "  0.381316    0.970574  -0.0144698       3.51786    3.3736     -24.4942  \n",
       "  1.75132     0.948956  -0.0165418       1.61441    4.22242    -22.6454  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = data[2:end,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## remaining columns characterize the density\n",
    "#basis_xgrid    = basis_logspline(xgrid,PhatKnots)\n",
    "#PDens_overtime   = zeros(size(data)[1],length(xgrid));\n",
    "#for tt = 1:size(data)[1]\n",
    " #   PDensCoef_tt    = coefRecover(data[tt,4:end]', \n",
    "   #                      PhatDensCoef_lambda, PhatDensCoef_mean)\n",
    "  # # PDensCoef_tt    = pdfNormalize_unrate(PDensCoef_tt+mean(seasonal_mean,1),PhatKnots,(data[tt,3] + mean(UNR)),0,3)       \n",
    "   # PDensCoef_tt    = pdfNormalize_unrate(PDensCoef_tt+mean(seasonal_mean,1),PhatKnots,0,0,3)  # use unemployment rate to get normalization constant\n",
    "   # PDens_overtime[tt,:] = exp.(basis_xgrid*PDensCoef_tt')'\n",
    "#end\n",
    "\n",
    "# remaining columns characterize the density\n",
    "basis_xgrid    = basis_logspline(xgrid,PhatKnots)\n",
    "PDens_overtime   = zeros(size(data)[1],length(xgrid));\n",
    "for tt = 1:size(data)[1]\n",
    "    PDensCoef_tt    = coefRecover(data[tt,4:end]', \n",
    "                         PhatDensCoef_lambda, PhatDensCoef_mean)\n",
    "    PDensCoef_tt    = pdfNormalize_unrate(PDensCoef_tt+mean(seasonal_mean,1),PhatKnots,(data[tt,3] + mean(UNR)),0,3)       \n",
    "    PDens_overtime[tt,:] = exp.(basis_xgrid*PDensCoef_tt')'\n",
    "    PDens_overtime[tt,1] = PDens_overtime[tt,1]+(data[tt,3]+mean(UNR))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0518662149122807"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean(UNR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#writecsv(savedir * \"PDens_overtime_mixedDist_seasonality.csv\", [period[2:end] PDens_overtime])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# remaining columns characterize the density\n",
    "basis_xgrid    = basis_logspline(xgrid,PhatKnots)\n",
    "PhatDens_overtime   = zeros(size(data)[1],length(xgrid));\n",
    "for tt = 1:size(data)[1]\n",
    "    PhatDensCoef_tt    = coefRecover(statepmean[tt,4:end]', \n",
    "                         PhatDensCoef_lambda, PhatDensCoef_mean)\n",
    "    PhatDensCoef_tt    = pdfNormalize_unrate(PhatDensCoef_tt+mean(seasonal_mean,1),PhatKnots,(statepmean[tt,3] + mean(UNR)),0,3)       \n",
    " #   PhatDensCoef_tt    = pdfNormalize_unrate(PhatDensCoef_tt+mean(seasonal_mean,1),PhatKnots,0,0,3)  # use unemployment rate to get normalization constant\n",
    "    PhatDens_overtime[tt,:] = exp.(basis_xgrid*PhatDensCoef_tt')'\n",
    "    PhatDens_overtime[tt,1] = PhatDens_overtime[tt,1]+(statepmean[tt,3]+mean(UNR))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#writecsv(savedir * \"PhatDens_overtime_mixedDist_seasonality.csv\", [period[2:end] PhatDens_overtime])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "function phatEval(x,statepmean_t,PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean)\n",
    "\n",
    "    PhatDensCoef_tt    = coefRecover(statepmean_t[4:end]', \n",
    "                         PhatDensCoef_lambda, PhatDensCoef_mean)\n",
    " #   PhatDensCoef_tt    = pdfNormalize_unrate(PhatDensCoef_tt+mean(seasonal_mean,1),PhatKnots,0,0,3)\n",
    "    PhatDensCoef_tt    = pdfNormalize_unrate(PhatDensCoef_tt+mean(seasonal_mean,1),PhatKnots, (statepmean_t[3] + mean(UNR)),0,3)\n",
    "                         #(statepmean_t[3] + mean(UNR)))  # use unemployment rate to get normalization constant \n",
    "    out                = exp.(basis_logspline(x,PhatKnots)*PhatDensCoef_tt')' # undo the transformation\n",
    "    return out\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "function phatIntegrate(statepmean_t,PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean, lb, ub)\n",
    "\n",
    "(Intout,error) = quadgk(x -> phatEval(x,statepmean_t,PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean),lb,ub)\n",
    "    return Intout\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "function phatIntegrate2(statepmean_t,PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean, lb, ub)\n",
    "\n",
    "(Intout,error) = quadgk(x -> phatEval(x,statepmean_t,PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean),lb,ub)\n",
    "    Intout = Intout + (statepmean_t[3]+mean(UNR))\n",
    "    return Intout\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "phatIntegrate(statepmean[1,:],PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean, 0.0, 0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "phatIntegrate2(statepmean[1,:],PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean, 0.0, 0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#statepmean_t = data[1,:]\n",
    "#phatIntegrate(statepmean_t, PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean, 0.0, 1.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ngrid = 20\n",
    "#grid_temp = linspace(0.2, log(1.8+1.0-0.2), ngrid)\n",
    "#grid_temp = exp.(grid_temp)-1\n",
    "grid_temp = linspace(0.2, 1.8, ngrid)\n",
    "grid_temp = [0.0; grid_temp]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "valuesIntegrate = zeros(length(grid_temp)-1,size(data)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## do the numerical integration to get CDF\n",
    "\n",
    "for tt = 1:size(data)[1]\n",
    "    statepmean_t = statepmean[tt,:]  # smoothed density\n",
    "\n",
    "#    statepmean_t = data[tt,:]  # raw density\n",
    "    for nn = 2:length(grid_temp)\n",
    "        if nn  == 2\n",
    "         valuesIntegrate[nn-1,tt] = phatIntegrate2(statepmean_t, PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean, grid_temp[nn-1], grid_temp[nn])[1]   \n",
    "        else\n",
    "         valuesIntegrate[nn-1,tt] = phatIntegrate(statepmean_t, PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean, grid_temp[nn-1], grid_temp[nn])[1]   \n",
    "        end            \n",
    "     print(\"grid # = $nn \\n\")\n",
    "    end\n",
    "     print(\"*************period = $tt \\n\")  \n",
    "end\n",
    "\n",
    "emp_cdf = [grid_temp[2:end] cumsum(valuesIntegrate,1)]\n",
    "\n",
    "#writecsv(savedir * \"emp_cdf_mixedDist2_seasonality.csv\", emp_cdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "writecsv(savedir * \"emp_cdf_mixedDist2_seasonality.csv\", emp_cdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## do the numerical integration to get CDF\n",
    "\n",
    "#for tt = 1:size(data)[1]\n",
    "#    statepmean_t = statepmean[tt,:]  # smoothed density\n",
    "##    statepmean_t = data[tt,:]  # raw density\n",
    "\n",
    "   # for nn = 2:length(grid_temp)\n",
    "   #  valuesIntegrate[nn-1,tt] = phatIntegrate(statepmean_t, PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, seasonal_mean, grid_temp[nn-1], grid_temp[nn])[1]   \n",
    "   #  print(\"grid # = $nn \\n\")\n",
    "   # end\n",
    "   #  print(\"*************period = $tt \\n\")  \n",
    "#end\n",
    "\n",
    "#emp_cdf = [grid_temp[2:end] cumsum(valuesIntegrate,1)]\n",
    "\n",
    "#writecsv(savedir * \"emp_cdf_mixedDist_seasonality.csv\", emp_cdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import Interpolations\n",
    "using Interpolations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now compute percentiles also taking into account of employment rate\n",
    "#unrate_t = data[:,3]+mean(UNR);\n",
    "unrate_t = statepmean[:,3]+mean(UNR);\n",
    "\n",
    "#statepmean = data\n",
    "\n",
    "emp_perc10 = zeros(size(statepmean)[1],1)\n",
    "emp_perc20 = zeros(size(statepmean)[1],1)\n",
    "emp_perc50 = zeros(size(statepmean)[1],1)\n",
    "emp_perc80 = zeros(size(statepmean)[1],1)\n",
    "emp_perc90 = zeros(size(statepmean)[1],1)\n",
    "\n",
    "for tt = 1:size(statepmean)[1]\n",
    "#tt = 1\n",
    " itp = interpolate((emp_cdf[:,tt+1],), emp_cdf[:,1], Gridded(Linear()))\n",
    "    emp_perc10[tt,1] = itp[0.1-unrate_t[tt]]\n",
    "    emp_perc20[tt,1] = itp[0.2-unrate_t[tt]]\n",
    "    emp_perc50[tt,1] = itp[0.5-unrate_t[tt]]\n",
    "    emp_perc80[tt,1] = itp[0.8-unrate_t[tt]]\n",
    "    emp_perc90[tt,1] = itp[0.9-unrate_t[tt]]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# undo the transformation\n",
    "\n",
    "emp_perc10 = 1/2*(exp.(emp_perc10)-exp.(-emp_perc10))\n",
    "emp_perc20 = 1/2*(exp.(emp_perc20)-exp.(-emp_perc20))\n",
    "emp_perc50 = 1/2*(exp.(emp_perc50)-exp.(-emp_perc50))\n",
    "emp_perc80 = 1/2*(exp.(emp_perc80)-exp.(-emp_perc80))\n",
    "emp_perc90 = 1/2*(exp.(emp_perc90)-exp.(-emp_perc90))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "emp_perc_unrate_added = [emp_perc10 emp_perc20 emp_perc50 emp_perc80 emp_perc90]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#writecsv(savedir * \"percentiles_unrate_added_raw_lbound_int1_seasonality.csv\", emp_perc_unrate_added)\n",
    "writecsv(savedir * \"percentiles_unrate_added_mixedDist2_seasonality.csv\", emp_perc_unrate_added)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "emp_perc10 = zeros(size(statepmean)[1],1)\n",
    "emp_perc20 = zeros(size(statepmean)[1],1)\n",
    "emp_perc50 = zeros(size(statepmean)[1],1)\n",
    "emp_perc80 = zeros(size(statepmean)[1],1)\n",
    "emp_perc90 = zeros(size(statepmean)[1],1)\n",
    "\n",
    "for tt = 1:size(statepmean)[1]\n",
    "#tt = 1\n",
    " itp = interpolate((emp_cdf[:,tt+1],), emp_cdf[:,1], Gridded(Linear()))\n",
    "    emp_perc10[tt,1] = itp[0.1]\n",
    "    emp_perc20[tt,1] = itp[0.2]\n",
    "    emp_perc50[tt,1] = itp[0.5]\n",
    "    emp_perc80[tt,1] = itp[0.8]\n",
    "    emp_perc90[tt,1] = itp[0.9]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# undo the transformation\n",
    "\n",
    "emp_perc10 = 1/2*(exp.(emp_perc10)-exp.(-emp_perc10))\n",
    "emp_perc20 = 1/2*(exp.(emp_perc20)-exp.(-emp_perc20))\n",
    "emp_perc50 = 1/2*(exp.(emp_perc50)-exp.(-emp_perc50))\n",
    "emp_perc80 = 1/2*(exp.(emp_perc80)-exp.(-emp_perc80))\n",
    "emp_perc90 = 1/2*(exp.(emp_perc90)-exp.(-emp_perc90))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "emp_perc = [emp_perc10 emp_perc20 emp_perc50 emp_perc80 emp_perc90]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "writecsv(savedir * \"percentiles_mixedDist2_seasonality.csv\", emp_perc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "statepmean[:,3]+mean(UNR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## without seasonal adjustment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "include(\"LoadData.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = data[2:end,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# remaining columns characterize the density\n",
    "basis_xgrid    = basis_logspline(xgrid,PhatKnots)\n",
    "PDens_overtime   = zeros(size(data)[1],length(xgrid));\n",
    "for tt = 1:size(data)[1]\n",
    "    PDensCoef_tt    = coefRecover(data[tt,4:end]', \n",
    "                         PhatDensCoef_lambda, PhatDensCoef_mean)\n",
    "    PDensCoef_tt    = pdfNormalize_unrate(PDensCoef_tt,PhatKnots,0,0,3)  # use unemployment rate to get normalization constant\n",
    "#    PDensCoef_tt    = pdfNormalize_unrate(PDensCoef_tt,PhatKnots,(data[tt,3] + mean(UNR)),0,3)  # use unemployment rate to get normalization constant\n",
    "    PDens_overtime[tt,:] = exp.(basis_xgrid*PDensCoef_tt')'\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## check the shape of the distribution\n",
    "dash = 0.3 * Compose.cm\n",
    "plot( layer(x=xgrid,y=PDens_overtime[1,:], Geom.line, \n",
    "            Theme(default_color=colorant\"red\",line_style=[dash] )),\n",
    "   # layer(x=xgrid,y=PhatDens_IRF[2,:], Geom.line, \n",
    "    #        Theme(default_color=colorant\"blue\",line_style=[dash] )),\n",
    "      Guide.xlabel(\"Period\"), Guide.ylabel(\"Earnings\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "function phatEval(x,statepmean_t,PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR)\n",
    "\n",
    "    PhatDensCoef_tt    = coefRecover(statepmean_t[4:end]', \n",
    "                         PhatDensCoef_lambda, PhatDensCoef_mean)\n",
    "    PhatDensCoef_tt    = pdfNormalize_unrate(PhatDensCoef_tt,PhatKnots,0,0,3)\n",
    "#    PhatDensCoef_tt    = pdfNormalize_unrate(PhatDensCoef_tt,PhatKnots,(statepmean_t[3] + mean(UNR)),0,3)\n",
    "                         #(statepmean_t[3] + mean(UNR)))  # use unemployment rate to get normalization constant \n",
    "    out                = exp.(basis_logspline(x,PhatKnots)*PhatDensCoef_tt')' # undo the transformation\n",
    "\n",
    "    return out\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "function phatIntegrate(statepmean_t,PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, lb, ub)\n",
    "\n",
    "(Intout,error) = quadgk(x -> phatEval(x,statepmean_t,PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR),lb,ub)\n",
    "    return Intout\n",
    "end\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ngrid = 20\n",
    "#grid_temp = linspace(0.2, log(1.8+1.0-0.2), ngrid)\n",
    "#grid_temp = exp.(grid_temp)-1\n",
    "grid_temp = linspace(0.2, 1.8, ngrid)\n",
    "grid_temp = [0.0; grid_temp]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "valuesIntegrate = zeros(length(grid_temp)-1,size(data)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## do the numerical integration to get CDF\n",
    "\n",
    "for tt = 1:10#size(data)[1]\n",
    "    statepmean_t = data[tt,:]  # raw density\n",
    "    for nn = 2:length(grid_temp)\n",
    "     valuesIntegrate[nn-1,tt] = phatIntegrate(statepmean_t, PhatDensCoef_lambda, PhatDensCoef_mean, PhatKnots, UNR, grid_temp[nn-1], grid_temp[nn])[1]   \n",
    "     print(\"grid # = $nn \\n\")\n",
    "    end\n",
    "     print(\"*************period = $tt \\n\")  \n",
    "end\n",
    "\n",
    "emp_cdf = [grid_temp[2:end] cumsum(valuesIntegrate,1)]\n",
    "#writecsv(savedir * \"emp_cdf_transformed.csv\", emp_cdf)\n",
    "#writecsv(savedir * \"emp_cdf_transformed_raw.csv\", emp_cdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import Interpolations\n",
    "using Interpolations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now compute percentiles also taking into account of employment rate\n",
    "unrate_t = data[:,3]+mean(UNR);\n",
    "\n",
    "statepmean = data\n",
    "\n",
    "emp_perc10 = zeros(size(statepmean)[1],1)\n",
    "emp_perc20 = zeros(size(statepmean)[1],1)\n",
    "emp_perc50 = zeros(size(statepmean)[1],1)\n",
    "emp_perc80 = zeros(size(statepmean)[1],1)\n",
    "emp_perc90 = zeros(size(statepmean)[1],1)\n",
    "\n",
    "for tt = 1:size(statepmean)[1]\n",
    "#tt = 1\n",
    " itp = interpolate((emp_cdf[:,tt+1],), emp_cdf[:,1], Gridded(Linear()))\n",
    "    emp_perc10[tt,1] = itp[0.1-unrate_t[tt]]\n",
    "    emp_perc20[tt,1] = itp[0.2-unrate_t[tt]]\n",
    "    emp_perc50[tt,1] = itp[0.5-unrate_t[tt]]\n",
    "    emp_perc80[tt,1] = itp[0.8-unrate_t[tt]]\n",
    "    emp_perc90[tt,1] = itp[0.9-unrate_t[tt]]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# undo the transformation\n",
    "\n",
    "emp_perc10 = 1/2*(exp.(emp_perc10)-exp.(-emp_perc10))\n",
    "emp_perc20 = 1/2*(exp.(emp_perc20)-exp.(-emp_perc20))\n",
    "emp_perc50 = 1/2*(exp.(emp_perc50)-exp.(-emp_perc50))\n",
    "emp_perc80 = 1/2*(exp.(emp_perc80)-exp.(-emp_perc80))\n",
    "emp_perc90 = 1/2*(exp.(emp_perc90)-exp.(-emp_perc90))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "emp_perc_unrate_added = [emp_perc10 emp_perc20 emp_perc50 emp_perc80 emp_perc90]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "writecsv(savedir * \"test_nolbound_int1_noseasonality.csv\", emp_perc_unrate_added[1:10,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#writecsv(savedir * \"empirical_percentiles_mat_unrate_added_raw_int1_noseason_nolbound.csv\", emp_perc_unrate_added[2:end,:])\n",
    "#writecsv(savedir * \"empirical_percentiles_mat_unrate_added_raw.csv\", emp_perc_unrate_added)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#writecsv(savedir * \"PDens_overtime_int1_noseason_nolbound.csv\", [period[2:end] PDens_overtime[2:end,:]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
